home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / offset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-14  |  27.6 KB  |  773 lines
  1. /******************************************************************************
  2. * Offset.c - computes offset approximation to curves and surfaces.            *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8. #include "extra_fn.h"
  9.  
  10. #define MAX_OFFSET_IMPROVE_ITERS    20
  11. #define NORMAL_PERTURB            1e-4
  12. #define OFFSET_TRIM_EPS            1.25
  13.  
  14. /*****************************************************************************
  15. * DESCRIPTION:                                                               M
  16. * Given a curve and an offset amount OffsetDist, returns an approximation to M
  17. * the offset curve by offseting the control polygon in the normal direction. M
  18. *                                                                            *
  19. * PARAMETERS:                                                                M
  20. *   Crv:          To approximate its offset curve with distance OffsetDist.  M
  21. *   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
  22. *   BezInterp:    If TRUE, control points are interpolated when the curve is M
  23. *          reduced to a Bezier form. Otherwise, control points are    M
  24. *          translated OffsetDist amount only, under estimating the    M
  25. *          Offset.                              M
  26. *                                                                            *
  27. * RETURN VALUE:                                                              M
  28. *   CagdCrvStruct *:   An approximation to the offset curve.                 M
  29. *                                                                            *
  30. * KEYWORDS:                                                                  M
  31. *   SymbCrvOffset, offset                                                    M
  32. *****************************************************************************/
  33. CagdCrvStruct *SymbCrvOffset(CagdCrvStruct *Crv,
  34.                  CagdRType OffsetDist,
  35.                  CagdBType BezInterp)
  36. {
  37.     CagdBType
  38.     IsBezierCurve = FALSE,
  39.     IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  40.     int i, j,
  41.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType),
  42.     Order = Crv -> Order,
  43.     Length = Crv -> Length;
  44.     CagdBType
  45.     HasNewKV = FALSE;
  46.     CagdVecStruct *N;
  47.     CagdCrvStruct
  48.     *OffsetCrv = CagdCrvCopy(Crv);
  49.     CagdRType *Nodes, *NodePtr,
  50.     **Points = OffsetCrv -> Points,
  51.     *KV = NULL;
  52.  
  53.     switch (Crv -> GType) {
  54.     case CAGD_CBEZIER_TYPE:
  55.         HasNewKV = TRUE;
  56.         KV = BspKnotUniformOpen(Length, Order, NULL);
  57.         IsBezierCurve = TRUE;
  58.         break;
  59.     case CAGD_CBSPLINE_TYPE:
  60.         KV = Crv -> KnotVector;
  61.         IsBezierCurve = Crv -> Length == Crv -> Order;
  62.         break;
  63.     case CAGD_CPOWER_TYPE:
  64.         SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
  65.         return NULL;
  66.     default:
  67.         SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
  68.         return NULL;
  69.     }
  70.  
  71.     Nodes = BspKnotNodes(KV, Length + Order, Order);
  72.     NodePtr = Nodes;
  73.  
  74.     /* Interpolate the compute control points instead of under estimating */
  75.     /* This offset.                                  */
  76.     if (BezInterp && IsBezierCurve) {
  77.     CagdCrvStruct *TCrv;
  78.     if (IsRational) {
  79.         TCrv = CagdCoerceCrvTo(OffsetCrv, CAGD_MAKE_PT_TYPE(FALSE,
  80.                    CAGD_NUM_OF_PT_COORD(OffsetCrv -> PType)));
  81.  
  82.         CagdCrvFree(OffsetCrv);
  83.         OffsetCrv = TCrv;
  84.         Points = OffsetCrv -> Points;
  85.     }
  86.  
  87.     for (j = 0; j < Length; j++, NodePtr++) {
  88.         CagdRType
  89.         *R = CagdCrvEval(Crv, *NodePtr);
  90.  
  91.         if ((N = CagdCrvNormal(Crv, *NodePtr)) == NULL &&
  92.         (N = CagdCrvNormal(Crv, NodePtr == Nodes ?
  93.                          *NodePtr + NORMAL_PERTURB :
  94.                          *NodePtr - NORMAL_PERTURB)) == NULL) {
  95.         SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL);
  96.         return NULL;
  97.         }
  98.  
  99.         for (i = 1; i <= MaxCoord; i++)
  100.         Points[i][j] = R[i] / (IsRational ? R[0] : 1.0) +
  101.                    N -> Vec[i - 1] * OffsetDist;
  102.     }
  103.  
  104.     TCrv = CagdCrvCopy(OffsetCrv);
  105.     for (i = 1; i <= MaxCoord; i++)
  106.         BzrCrvInterp(TCrv -> Points[i], OffsetCrv -> Points[i], Length);
  107.  
  108.     CagdCrvFree(OffsetCrv);
  109.     OffsetCrv = TCrv;
  110.     }
  111.     else {
  112.     for (j = 0; j < Length; j++, NodePtr++) {
  113.         if ((N = CagdCrvNormal(Crv, *NodePtr)) == NULL &&
  114.         (N = CagdCrvNormal(Crv, NodePtr == Nodes ?
  115.                          *NodePtr + NORMAL_PERTURB :
  116.                          *NodePtr - NORMAL_PERTURB)) == NULL) {
  117.         SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL);
  118.         return NULL;
  119.         }
  120.  
  121.         for (i = 1; i <= MaxCoord; i++)
  122.         Points[i][j] +=
  123.             N -> Vec[i - 1] * OffsetDist *
  124.             (IsRational ? Points[W][j] : 1.0);
  125.     }
  126.     }
  127.  
  128.     if (HasNewKV)
  129.     IritFree((VoidPtr) KV);
  130.     IritFree((VoidPtr) Nodes);
  131.  
  132.     return OffsetCrv;
  133. }
  134.  
  135. /*****************************************************************************
  136. * DESCRIPTION:                                                               M
  137. * Given a curve and an offset amount OffsetDist, returns an approximation to M
  138. * the offset curve by offseting the control polygon in the normal direction. M
  139. *   If resulting offset is not satisfying the required tolerance the curve   M
  140. * is subdivided and the algorithm recurses on both parts.             M
  141. *                                                                            *
  142. * PARAMETERS:                                                                M
  143. *   Crv:          To approximate its offset curve with distance OffsetDist.  M
  144. *   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
  145. *   Tolerance:    Accuracy control.                                          M
  146. *   BezInterp:    If TRUE, control points are interpolated when the curve is M
  147. *          reduced to a Bezier form. Otherwise, control points are    M
  148. *          translated OffsetDist amount only, under estimating the    M
  149. *          Offset.                              M
  150. *                                                                            *
  151. * RETURN VALUE:                                                              M
  152. *   CagdCrvStruct *:   An approximation to the offset curve, to within       M
  153. *                      Tolerance.                         M
  154. *                                                                            *
  155. * KEYWORDS:                                                                  M
  156. *   SymbCrvSubdivOffset, offset                                              M
  157. *****************************************************************************/
  158. CagdCrvStruct *SymbCrvSubdivOffset(CagdCrvStruct *Crv,
  159.                    CagdRType OffsetDist,
  160.                    CagdRType Tolerance,
  161.                    CagdBType BezInterp)
  162. {
  163.     CagdCrvStruct
  164.     *OffCrv = SymbCrvOffset(Crv, OffsetDist, BezInterp),
  165.     *Dist = SymbCrvSub(Crv, OffCrv),
  166.     *DistSqr = SymbCrvDotProd(Dist, Dist);
  167.     CagdRType *R, MinVal, MaxVal, TMin, TMax;
  168.  
  169.     CagdCrvFree(Dist);
  170.  
  171.     R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, TRUE);
  172.     MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
  173.     R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, FALSE);
  174.     MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
  175.  
  176.     CagdCrvFree(DistSqr);
  177.  
  178.     CagdCrvDomain(Crv, &TMin, &TMax);
  179.     if ((ABS(MinVal - ABS(OffsetDist)) > Tolerance ||
  180.      ABS(MaxVal - ABS(OffsetDist)) > Tolerance) &&
  181.     TMax - TMin > NORMAL_PERTURB * 10.0) {
  182.     CagdBType NewCrv;
  183.     CagdCrvStruct *Crv1, *Crv2, *OffCrv1, *OffCrv2;
  184.  
  185.     if (CAGD_IS_BEZIER_CRV(Crv)) {
  186.         /* Promote to a Bspline so we can keep track of parametrization. */
  187.         Crv = CnvrtBezier2BsplineCrv(Crv);
  188.         NewCrv = TRUE;
  189.     }
  190.     else
  191.         NewCrv = FALSE;
  192.  
  193.     if (TMax - TMin > NORMAL_PERTURB * 10.0) {
  194.         CagdCrvFree(OffCrv);
  195.  
  196.         Crv1 = CagdCrvSubdivAtParam(Crv, (TMin + TMax) / 2.0);
  197.         Crv2 = Crv1 -> Pnext;
  198.         Crv1 -> Pnext = NULL;
  199.         OffCrv1 = SymbCrvSubdivOffset(Crv1, OffsetDist, Tolerance, BezInterp);
  200.         OffCrv2 = SymbCrvSubdivOffset(Crv2, OffsetDist, Tolerance, BezInterp);
  201.         CagdCrvFree(Crv1);
  202.         CagdCrvFree(Crv2);
  203.  
  204.         OffCrv = CagdMergeCrvCrv(OffCrv1, OffCrv2, TRUE);
  205.         CagdCrvFree(OffCrv1);
  206.         CagdCrvFree(OffCrv2);
  207.     }
  208.  
  209.     if (NewCrv)
  210.         CagdCrvFree(Crv);
  211.     }
  212.  
  213.     return OffCrv;
  214. }
  215.  
  216. /*****************************************************************************
  217. * DESCRIPTION:                                                               M
  218. * Given a surface and an offset amount OffsetDist, returns an approximation  M
  219. * to the offset surface by offseting the control mesh in the normal         M
  220. * direction.                                                                 M
  221. *                                                                            *
  222. * PARAMETERS:                                                                M
  223. *   Srf:         To approximate its offset surface with distance OffsetDist. M
  224. *   OffsetDist:  Amount of offset. Negative denotes other offset direction.  M
  225. *                                                                            *
  226. * RETURN VALUE:                                                              M
  227. *   CagdSrfStruct *:   An approximation to the offset surface.               M
  228. *                                                                            *
  229. * KEYWORDS:                                                                  M
  230. *   SymbSrfOffset, offset                                                    M
  231. *****************************************************************************/
  232. CagdSrfStruct *SymbSrfOffset(CagdSrfStruct *Srf, CagdRType OffsetDist)
  233. {
  234.     CagdBType
  235.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  236.     int    i, Row, Col,
  237.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType),
  238.     UOrder = Srf -> UOrder,
  239.     VOrder = Srf -> VOrder,
  240.     ULength = Srf -> ULength,
  241.     VLength = Srf -> VLength;
  242.     CagdBType
  243.     HasNewKV = FALSE;
  244.     CagdVecStruct *N;
  245.     CagdSrfStruct
  246.     *OffsetSrf = CagdSrfCopy(Srf);
  247.     CagdRType *UNodes, *VNodes, *UNodePtr, *VNodePtr,
  248.     **Points = OffsetSrf -> Points,
  249.     *UKV = NULL,
  250.     *VKV = NULL;
  251.  
  252.     switch (Srf -> GType) {
  253.     case CAGD_SBEZIER_TYPE:
  254.         HasNewKV = TRUE;
  255.         UKV = BspKnotUniformOpen(ULength, UOrder, NULL);
  256.         VKV = BspKnotUniformOpen(VLength, VOrder, NULL);
  257.         break;
  258.     case CAGD_SBSPLINE_TYPE:
  259.         UKV = Srf -> UKnotVector;
  260.         VKV = Srf -> VKnotVector;
  261.         break;
  262.     case CAGD_SPOWER_TYPE:
  263.         SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
  264.         return NULL;
  265.     default:
  266.         SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
  267.         return NULL;
  268.     }
  269.  
  270.     UNodes = BspKnotNodes(UKV, ULength + UOrder, UOrder);
  271.     VNodes = BspKnotNodes(VKV, VLength + VOrder, VOrder);
  272.  
  273.     if (IsNotRational)
  274.     for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++)
  275.         for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) {
  276.             N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr);
  277.             for (i = 1; i <= MaxCoord; i++)
  278.             Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] +=
  279.             N -> Vec[i - 1] * OffsetDist;
  280.         }
  281.     else
  282.     for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++)
  283.         for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) {
  284.             N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr);
  285.             for (i = 1; i <= MaxCoord; i++)
  286.                 Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] +=
  287.             N -> Vec[i - 1] * OffsetDist *
  288.                 Points[W][CAGD_MESH_UV(OffsetSrf, Col, Row)];
  289.         }
  290.  
  291.     if (HasNewKV) {
  292.     IritFree((VoidPtr) UKV);
  293.     IritFree((VoidPtr) VKV);
  294.     }
  295.  
  296.     IritFree((VoidPtr) UNodes);
  297.     IritFree((VoidPtr) VNodes);
  298.  
  299.     return OffsetSrf;
  300. }
  301.  
  302. /*****************************************************************************
  303. * DESCRIPTION:                                                               M
  304. * Given a surface and an offset amount OffsetDist, returns an approximation  M
  305. * to the offset surface by offseting the control mesh in the normal         M
  306. * direction.                                                                 M
  307. *   If resulting offset is not satisfying the required tolerance the surface M
  308. * is subdivided and the algorithm recurses on both parts.             M
  309. *                                                                            *
  310. * PARAMETERS:                                                                M
  311. *   Srf:         To approximate its offset surface with distance OffsetDist. M
  312. *   OffsetDist:  Amount of offset. Negative denotes other offset direction.  M
  313. *   Tolerance:    Accuracy control.                                          M
  314. *                                                                            *
  315. * RETURN VALUE:                                                              M
  316. *   CagdSrfStruct *:   An approximation to the offset surface, to within     M
  317. *                      Tolerance.                             M
  318. *                                                                            *
  319. * KEYWORDS:                                                                  M
  320. *   SymbSrfSubdivOffset, offset                                              M
  321. *****************************************************************************/
  322. CagdSrfStruct *SymbSrfSubdivOffset(CagdSrfStruct *Srf,
  323.                    CagdRType OffsetDist,
  324.                    CagdRType Tolerance)
  325. {
  326.     CagdSrfStruct
  327.     *OffSrf = SymbSrfOffset(Srf, OffsetDist),
  328.     *Dist = SymbSrfSub(Srf, OffSrf),
  329.     *DistSqr = SymbSrfDotProd(Dist, Dist);
  330.     CagdRType *R, MinVal, MaxVal;
  331.     CagdSrfDirType Dir;
  332.  
  333.     CagdSrfFree(Dist);
  334.  
  335.     R = SymbExtremumCntPtVals(DistSqr -> Points,
  336.                   DistSqr -> ULength * DistSqr -> VLength, TRUE);
  337.     MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
  338.     R = SymbExtremumCntPtVals(DistSqr -> Points,
  339.                   DistSqr -> ULength * DistSqr -> VLength, FALSE);
  340.     MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
  341.  
  342.     CagdSrfFree(DistSqr);
  343.  
  344.     if (ABS(MinVal - ABS(OffsetDist)) > Tolerance ||
  345.     ABS(MaxVal - ABS(OffsetDist)) > Tolerance) {
  346.     CagdSrfStruct *Srf1, *Srf2, *OffSrf1, *OffSrf2;
  347.     CagdRType UMin, UMax, VMin, VMax, t;
  348.     CagdBType NewSrf;
  349.  
  350.     /* Make sure it is a Bspline surface. */
  351.     if (CAGD_IS_BEZIER_SRF(Srf)) {
  352.         /* Promote to a Bspline so we can keep track of parametrization. */
  353.         Srf = CnvrtBezier2BsplineSrf(Srf);
  354.         NewSrf = TRUE;
  355.     }
  356.     else
  357.         NewSrf = FALSE;
  358.  
  359.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  360.     if (MAX(UMax - UMin, VMax - VMin) > NORMAL_PERTURB * 10.0) {
  361.         CagdSrfFree(OffSrf);
  362.  
  363.         if (UMax - UMin > VMax - VMin) {
  364.         Dir = CAGD_CONST_U_DIR;
  365.         t = (UMin + UMax) / 2.0;
  366.         }
  367.         else {
  368.         Dir = CAGD_CONST_V_DIR;
  369.         t = (VMin + VMax) / 2.0;
  370.         }
  371.         Srf1 = CagdSrfSubdivAtParam(Srf, t, Dir);
  372.         Srf2 = Srf1 -> Pnext;
  373.         Srf1 -> Pnext = NULL;
  374.         OffSrf1 = SymbSrfSubdivOffset(Srf1, OffsetDist, Tolerance);
  375.         OffSrf2 = SymbSrfSubdivOffset(Srf2, OffsetDist, Tolerance);
  376.         CagdSrfFree(Srf1);
  377.         CagdSrfFree(Srf2);
  378.  
  379.         OffSrf = CagdMergeSrfSrf(OffSrf1, OffSrf2, Dir, TRUE, TRUE);
  380.         CagdSrfFree(OffSrf1);
  381.         CagdSrfFree(OffSrf2);
  382.     }
  383.  
  384.     if (NewSrf)
  385.         CagdSrfFree(Srf);
  386.     }
  387.  
  388.     return OffSrf;
  389. }
  390.  
  391. /*****************************************************************************
  392. * DESCRIPTION:                                                               M
  393. * Given a curve and an offset amount OffsetDist, returns an approximation to M
  394. * the offset curve by offseting the control polygon in the normal direction. M
  395. *   This function computes an approximation to the offset using             M
  396. * OffsetAprxFunc, measure the error and use it to refine and decrease the    M
  397. * error adaptively.                                 M
  398. *   Bezier curves are promoted to Bsplines curves.                 M
  399. *   See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable         M
  400. * Distance Offset Operator for Free Form Curves and Surfaces". International M
  401. * Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March    M
  402. * 1991, pp 67-78.                                 M
  403. *                                                                            *
  404. * PARAMETERS:                                                                M
  405. *   OrigCrv:      To approximate its offset curve with distance OffsetDist.  M
  406. *   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
  407. *   OffsetError:  Tolerance control.                                         M
  408. *   OffsetAprxFunc:  A function that can be used to approximate an offset    M
  409. *                    of a curve. If NULL SymbCrvOffset function is selected. M
  410. *   BezInterp:    If TRUE, control points are interpolated when the curve is M
  411. *          reduced to a Bezier form. Otherwise, control points are    M
  412. *          translated OffsetDist amount only, under estimating the    M
  413. *          Offset.                              M
  414. *                                                                            *
  415. * RETURN VALUE:                                                              M
  416. *   CagdCrvStruct *:   An approximation to the offset curve, to within       M
  417. *                      OffsetError.                                          M
  418. *                                                                            *
  419. * KEYWORDS:                                                                  M
  420. *   SymbCrvAdapOffset, offset                                                M
  421. *****************************************************************************/
  422. CagdCrvStruct *SymbCrvAdapOffset(CagdCrvStruct *OrigCrv,
  423.                  CagdRType OffsetDist,
  424.                  CagdRType OffsetError,
  425.                  SymbOffCrvFuncType OffsetAprxFunc,
  426.                  CagdBType BezInterp)
  427. {
  428.     int i;
  429.     CagdBType
  430.     IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv);
  431.     CagdRType Min, Max, TMin, TMax,
  432.     OffsetDistSqr = SQR(OffsetDist);
  433.     CagdCrvStruct *OffsetCrv, *Crv;
  434.  
  435.     switch (OrigCrv -> GType) {
  436.     case CAGD_CBEZIER_TYPE:
  437.         Crv = CnvrtBezier2BsplineCrv(OrigCrv);
  438.         break;
  439.     case CAGD_CBSPLINE_TYPE:
  440.         Crv = CagdCrvCopy(OrigCrv);
  441.         break;
  442.     case CAGD_CPOWER_TYPE:
  443.     default:
  444.         SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
  445.         Crv = NULL;
  446.         break;
  447.     }
  448.  
  449.     if (OffsetAprxFunc == NULL)
  450.     OffsetAprxFunc = SymbCrvOffset;
  451.  
  452.     CagdCrvDomain(Crv, &TMin, &TMax);
  453.  
  454.     for (i = 0; i < MAX_OFFSET_IMPROVE_ITERS; i++) {
  455.     CagdCrvStruct *DiffCrv, *DistSqrCrv;
  456.  
  457.     OffsetCrv = OffsetAprxFunc(Crv, OffsetDist, BezInterp);
  458.     DiffCrv = SymbCrvSub(OffsetCrv, Crv);
  459.     DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv);
  460.     CagdCrvFree(DiffCrv);
  461.  
  462.     CagdCrvMinMax(DistSqrCrv, 1, &Min, &Max);
  463.  
  464.     if (OffsetDistSqr - Min < OffsetError &&
  465.         Max - OffsetDistSqr < OffsetError) {
  466.         /* Error is within bounds - returns this offset approximation. */
  467.         CagdCrvFree(DistSqrCrv);
  468.         break;
  469.     }
  470.     else {
  471.         /* Refine in regions where the error is too large. */
  472.         int j, k,
  473.             Length = DistSqrCrv -> Length,
  474.             Order = DistSqrCrv -> Order,
  475.             KVLen = Length + Order;
  476.         CagdRType
  477.         *KV = DistSqrCrv -> KnotVector,
  478.         *Nodes = BspKnotNodes(KV, KVLen, Order),
  479.         *RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * 2 * Length);
  480.  
  481.         for (j = k = 0; j < Length; j++) {
  482.         CagdRType
  483.             *Pt = CagdCrvEval(DistSqrCrv, Nodes[j]),
  484.             V = OffsetDistSqr - (IsRational ? Pt[1] / Pt[0] : Pt[1]);
  485.  
  486.         if (ABS(V) > OffsetError) {
  487.             int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);
  488.  
  489.             if (APX_EQ(KV[Index], Nodes[j])) {
  490.             if (j > 0)
  491.                 RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
  492.             if (j < Length - 1)
  493.                 RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
  494.             }
  495.             else
  496.             RefKV[k++] = Nodes[j];
  497.         }
  498.         }
  499.         CagdCrvFree(DistSqrCrv);
  500.         IritFree((VoidPtr) Nodes);
  501.  
  502.         if (k == 0) {
  503.         /* No refinement was found needed - return current curve. */
  504.         IritFree((VoidPtr) RefKV);
  505.         break;
  506.         }
  507.         else {
  508.         CagdCrvStruct
  509.             *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);
  510.  
  511.         IritFree((VoidPtr) RefKV);
  512.         CagdCrvFree(Crv);
  513.         Crv = CTmp;
  514.         }
  515.     }
  516.     }
  517.  
  518.     CagdCrvFree(Crv);
  519.     return OffsetCrv;
  520. }
  521.  
  522. /*****************************************************************************
  523. * DESCRIPTION:                                                               M
  524. * Same function as CagdCrvAdapOffset, but trims the self intersection loops. M
  525. *   See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable         M
  526. * Distance Offset Operator for Free Form Curves and Surfaces". International M
  527. * Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March    M
  528. * 1991, pp 67-78.                                 M
  529. *                                                                            *
  530. * PARAMETERS:                                                                M
  531. *   OrigCrv:      To approximate its offset curve with distance OffsetDist.  M
  532. *   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
  533. *   OffsetError:  Tolerance control.                                         M
  534. *   OffsetAprxFunc:  A function that can be used to approximate an offset    M
  535. *                    of a curve. If NULL SymbCrvOffset function is selected. M
  536. *             Third parameter of SymbOffCrvFuncType is optional.         M
  537. *   BezInterp:    If TRUE, control points are interpolated when the curve is M
  538. *          reduced to a Bezier form. Otherwise, control points are    M
  539. *          translated OffsetDist amount only, under estimating the    M
  540. *          Offset.                              M
  541. *                                                                            *
  542. * RETURN VALUE:                                                              M
  543. *   CagdCrvStruct *:   An approximation to the offset curve, to within       M
  544. *                      OffsetError.                                          M
  545. *                                                                            *
  546. * KEYWORDS:                                                                  M
  547. *   SymbCrvAdapOffsetTrim, offset                                            M
  548. *****************************************************************************/
  549. CagdCrvStruct *SymbCrvAdapOffsetTrim(CagdCrvStruct *OrigCrv,
  550.                      CagdRType OffsetDist,
  551.                      CagdRType OffsetError,
  552.                      SymbOffCrvFuncType OffsetAprxFunc,
  553.                      CagdBType BezInterp)
  554. {
  555.     CagdBType
  556.         IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv);
  557.     CagdCrvStruct *CrvtrSqr, *CrvtrSign, *Crv, *Crvs, *LastCrv,
  558.     *CrvOffsetList = NULL;
  559.     CagdPtStruct *Pts, *PtsHead, *LastPts;
  560.     CagdRType
  561.         OffsetDistSqr1 = 1.0 / SQR(OffsetDist);
  562.  
  563.     switch (OrigCrv -> GType) {
  564.     case CAGD_CBEZIER_TYPE:
  565.         Crv = CnvrtBezier2BsplineCrv(OrigCrv);
  566.         break;
  567.     case CAGD_CBSPLINE_TYPE:
  568.         Crv = CagdCrvCopy(OrigCrv);
  569.         break;
  570.     default:
  571.         SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
  572.         Crv = NULL;
  573.         break;
  574.     }
  575.  
  576.     if (OffsetDist == 0.0)
  577.     return Crv;
  578.  
  579.     CrvtrSqr = SymbCrv2DCurvatureSqr(Crv);
  580.     CrvtrSign = SymbCrv2DCurvatureSign(Crv);
  581.  
  582.     PtsHead = SymbCrvConstSet(CrvtrSqr, 1, SQR(OffsetError),
  583.                   OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS));
  584.  
  585.     for (Pts = PtsHead, LastPts = NULL; Pts != NULL;) {
  586.     CagdRType *CrvtrSignPt, CrvtrSignVal;
  587.  
  588.     CrvtrSignPt = CagdCrvEval(CrvtrSign, Pts -> Pt[0]);
  589.     CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0]
  590.                   : CrvtrSignPt[1];
  591.  
  592.     if (CrvtrSignVal * OffsetDist > 0.0) {
  593.         /* Remove point from list. */
  594.         if (LastPts != NULL) {
  595.         LastPts -> Pnext = Pts -> Pnext;
  596.         CagdPtFree(Pts);
  597.         Pts = LastPts -> Pnext;
  598.         }
  599.         else {
  600.         PtsHead = PtsHead -> Pnext;
  601.         CagdPtFree(Pts);
  602.         Pts = PtsHead;
  603.         }
  604.     }
  605.     else {
  606.         LastPts = Pts;
  607.         Pts = Pts -> Pnext;
  608.     }
  609.     }
  610.  
  611.     for (Pts = PtsHead;
  612.      Pts != NULL || Crv != NULL;
  613.      Pts = (Pts != NULL ? Pts -> Pnext : NULL)) {
  614.     CagdCrvStruct
  615.         *Crv1 = Pts ? CagdCrvSubdivAtParam(Crv, Pts -> Pt[0])
  616.             : CagdCrvCopy(Crv),
  617.         *Crv2 = Pts ? Crv1 -> Pnext : NULL;
  618.     CagdRType Min, Max, *CrvtrSqrPt, CrvtrSqrVal,
  619.                         *CrvtrSignPt, CrvtrSignVal;
  620.  
  621.     CagdCrvDomain(Crv1, &Min, &Max);
  622.  
  623.     CrvtrSqrPt = CagdCrvEval(CrvtrSqr, (Min + Max) / 2);
  624.     CrvtrSqrVal = CrvtrSqrPt[1] / CrvtrSqrPt[0];
  625.     CrvtrSignPt = CagdCrvEval(CrvtrSign, (Min + Max) / 2);
  626.     CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0]
  627.                   : CrvtrSignPt[1];
  628.  
  629.     if (CrvtrSqrVal < OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS) ||
  630.         CrvtrSignVal * OffsetDist > 0.0) {
  631.         CagdCrvStruct
  632.         *Crv1Off = SymbCrvAdapOffset(Crv1, OffsetDist, OffsetError,
  633.                          OffsetAprxFunc, BezInterp);
  634.  
  635.         LIST_PUSH(Crv1Off, CrvOffsetList);
  636.     }
  637.  
  638.     CagdCrvFree(Crv1);
  639.     CagdCrvFree(Crv);
  640.     Crv = Crv2;
  641.     }
  642.  
  643.     Crvs = CagdListReverse(CrvOffsetList);
  644.     CrvOffsetList = NULL;
  645.     LastCrv = Crvs;
  646.     Crvs = Crvs -> Pnext;
  647.     for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) {
  648.     CagdCrvStruct *CTmp, *CTmp2;
  649.     CagdSrfStruct
  650.         *DistSrf = SymbSrfDistCrvCrv(LastCrv, Crv);
  651.     CagdPtStruct
  652.         *IPts = SymbSrfDistFindPoints(DistSrf, OffsetError, FALSE);
  653.  
  654.     CagdSrfFree(DistSrf);
  655.  
  656.     if (IPts != NULL) {
  657.         if (IPts -> Pnext) {
  658.         SYMB_FATAL_ERROR(SYMB_ERR_TOO_COMPLEX);
  659.         return NULL;
  660.         }
  661.         else {
  662.         CagdRType TMin, TMax;
  663.  
  664.         CagdCrvDomain(LastCrv, &TMin, &TMax);
  665.         CTmp = CagdCrvRegionFromCrv(LastCrv, TMin, IPts -> Pt[0]);
  666.         CagdCrvFree(LastCrv);
  667.         LastCrv = CTmp;
  668.  
  669.         CagdCrvDomain(Crv, &TMin, &TMax);
  670.         CTmp = CagdCrvRegionFromCrv(Crv, IPts -> Pt[1], TMax);
  671.  
  672.         CTmp2 = CagdMergeCrvCrv(LastCrv, CTmp, FALSE);
  673.         CagdCrvFree(CTmp);
  674.         CagdCrvFree(LastCrv);
  675.         LastCrv = CTmp2;
  676.         }
  677.     }
  678.     else {
  679.         /* Simply chain the pieces together. */
  680.         CTmp = CagdMergeCrvCrv(LastCrv, Crv, FALSE);
  681.         CagdCrvFree(LastCrv);
  682.         LastCrv = CTmp;
  683.     }
  684.     }
  685.  
  686.     CagdPtFreeList(PtsHead);
  687.     CagdCrvFree(CrvtrSqr);
  688.     CagdCrvFree(CrvtrSign);
  689.  
  690.     return LastCrv;
  691. }
  692.  
  693. /*****************************************************************************
  694. * DESCRIPTION:                                                               M
  695. * Given a curve and an offset amount OffsetDist, returns an approximation to M
  696. * the offset curve by least square fitting a curve to samples taken on the   M
  697. * offset curve.                                     M
  698. *   Resulting curve of order Order (degree of Crv if Order == 0) will have   M
  699. * NumOfDOF control points that least sqaure fit NumOfSamples samples on the  M
  700. * offset curve.                                     M
  701. *   Tolerance will be updated to hold an error distance measure.         M
  702. *                                                                            *
  703. * PARAMETERS:                                                                M
  704. *   Crv:          To approximate its offset curve with distance OffsetDist.  M
  705. *   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
  706. *   NumOfSamples: Number of samples to sample the offset curve at.           M
  707. *   NumOfDOF:     Number of degrees of freedom on the newly computed offset  M
  708. *          approximation. This is thesame as the number of control    M
  709. *          points the new curve will have.                 M
  710. *   Order:        Of the newly constructed offset approximation. If equal to M
  711. *          zero, the order of Crv will be used.                 M
  712. *   Tolerance:    To return an error estimate in the L-infinity norm.        M
  713. *                                                                            *
  714. * RETURN VALUE:                                                              M
  715. *   CagdCrvStruct *:   An approximation to the offset curve.             M
  716. *                                                                            *
  717. * KEYWORDS:                                                                  M
  718. *   SymbCrvLeastSquarOffset, offset                                          M
  719. *****************************************************************************/
  720. CagdCrvStruct *SymbCrvLeastSquarOffset(CagdCrvStruct *Crv,
  721.                        CagdRType OffsetDist,
  722.                        int NumOfSamples,
  723.                        int NumOfDOF,
  724.                        int Order,
  725.                        CagdRType *Tolerance)
  726. {
  727.     int i;
  728.     CagdRType TMin, TMax, t, dt;
  729.     CagdVType Tangent;
  730.     CagdCrvStruct *OffApprox,
  731.     *TangentCrv = CagdCrvDerive(Crv);
  732.     CagdPtStruct
  733.     *Pt = NULL,
  734.     *PtList = NULL;
  735.  
  736.     /* Create NumOfSamples samples on the offset curve. */
  737.     CagdCrvDomain(Crv, &TMin, &TMax);
  738.     dt = (TMax - TMin) / (NumOfSamples - 1);
  739.     for (i = 0, t = TMin; i < NumOfSamples; i++, t += dt) {
  740.     CagdRType *R;
  741.  
  742.     if (PtList == NULL)
  743.         PtList = Pt = CagdPtNew();
  744.     else {
  745.         Pt -> Pnext = CagdPtNew();
  746.         Pt = Pt -> Pnext;
  747.     }
  748.  
  749.     R = CagdCrvEval(Crv, t);
  750.     CagdCoerceToE3(Pt -> Pt, &R, -1, Crv -> PType);
  751.  
  752.     R = CagdCrvEval(TangentCrv, t);
  753.     CagdCoerceToE2(Tangent, &R, -1, TangentCrv -> PType);
  754.     Tangent[2] = 0.0;
  755.     PT_NORMALIZE(Tangent);
  756.     
  757.     Pt -> Pt[0] += Tangent[1] * OffsetDist;
  758.     Pt -> Pt[1] -= Tangent[0] * OffsetDist;
  759.     }
  760.  
  761.     OffApprox = BspCrvInterpPts(PtList,
  762.                 Order == 0 ? Crv -> Order : Order,
  763.                 MIN(NumOfDOF, NumOfSamples),
  764.                 CAGD_UNIFORM_PARAM);
  765.  
  766.     *Tolerance = BspCrvInterpPtsError(OffApprox, PtList, CAGD_UNIFORM_PARAM);
  767.  
  768.     CagdPtFreeList(PtList);
  769.     CagdCrvFree(TangentCrv);
  770.  
  771.     return OffApprox;
  772. }
  773.